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Abstract We pose and solve the analogue of Slepian's time-frequency concentration 
problem in the two-dimensional plane, for applications in the natural sciences. We 
determine an orthogonal family of strictly bandlimited functions that are optimally 
concentrated within a closed region of the plane, or, alternatively, of strictly space- 
limited functions that are optimally concentrated in the Fourier domain. The Carte- 
413 sian Slepian functions can be found by solving a Fredholm integral equation whose 

^ associated eigenvalues are a measure of the spatiospectral concentration. Both the 

Ch spatial and spectral regions of concentration can, in principle, have arbitrary geom- 

■ etry. However, for practical applications of signal representation or spectral analysis 

such as exist in geophysics or astronomy, in physical space irregular shapes, and in 
spectral space symmetric domains will usually be preferred. When the concentra- 
tion domains are circularly symmetric in both spaces, the Slepian functions are also 
^^ eigenfunctions of a Sturm-Liouville operator, leading to special algorithms for this 

fS| case, as is well known. Much like their one-dimensional and spherical counterparts 

^/^ with which we discuss them in a common framework, a basis of functions that are 

f — . simultaneously spatially and spectrally localized on arbitrary Cartesian domains will 

^^ be of great utility in many scientific disciplines, but especially in the geosciences. 
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1 Introduction 



The one-dimensional prolate spheroidal wave functions (pswf) have enjoyed an en- 
during popularity in the signal processing community ever since their introduction in 
the early 1960s ( Landau and PollaJ^l 961^1 962^ ^Slepian and Pollak^ 1961] ). Indeed, in 
many scientific and engineering disciplines the pswf and their relatives the discrete 
prolate spheroidal sequences (dpss) \GrUnbauw}\l9'^\\ \Slepian^91'^) have by now 
become the preferred data windows to regularize the quadratic inverse problem of 



power spectral estimation from time-series observations of finite extent (Percival and 



Waldeni l993). At the deliberate cost of introducing spectral bias, windowing the data 
with a set of such orthogonal "tapers" lowers the variance of the "multitaper" aver- 
age ^Thomson\19^2^ , which results in estimates of the power spectral density that are 
low-error in the mean-squared sense ^Thomson\\1990\. As a basis for function rep- 
resentation, approximation and interpolation jDe lsarte et al.\19^5\ Moore and Cada 
[2004; Shkolnisk y et al.\2{){)6\\Xiao et al 2001 ), or in stochastic linear inverse prob- 
lems \Bertero et <2/.|1985a|b[|Jg Villiers et al\2S^\WinghamH 992), the pswf have 



been less in the public eye, especially compared to wavelet analysis \Daubechies 



1992[ \Percival and Waldew\20Q6) , though, due to advances in computation, there 
has been a resurgent interest in recent years (Beylkin and Monzon 2002; \Karoui and 



Moumni\200^\\Khare and George\2003\\Walter and Sole ski ^2{){) 5), in particular as re- 



lates to using them for the numerical solution of partial differential equations \Beylkin 
and Sandberg\2m5\^oyd^m^^M^\Chen et k|2005| ) 



The pswf are the solutions to what has come to be known as "the concentration 
problem" ^Flandrm\l99^\\Percival and Walden\\99?>) of Slepian, Landau and Pollak, 
in which the energy of a bandlimited function is maximized by quadratic optimization 
inside a certain interval of time. Vice versa, it refers to the maximization of the spec- 
tral localization of a timelimited function inside a certain target bandwidth {Slepian 



1983| ). In the first version of the problem, the bandlimited, time-concentrated pswf 



form an orthogonal basis for the entire space of bandlimited signals that is also or- 
thogonal over the particular time interval of interest. In the second version the time- 
limited, band-concentrated pswf are a basis for square-integrable broadband signals 
that are exactly confined to the interval { Landau and Pollak\l96l\ Slepian and Pol 
7a^|1961| ). In general we shall refer to all singular functions of time-bandwidth or 



space-bandwidth projection operators as "Slepian functions". 

The fixed prescription of the "region of interest" in physical or spectral space 
is a deliberately narrow point of view that is well suited to scientific or engineer- 
ing studies where the assumption of stationarity, prior information, or the availabil- 
ity of data will dictate the interval of study from the outset. This distinguishes the 
Slepian functions philosophically from the eigenf unctions of full-phase- space local- 



ization operators {Daubechies 1988 ; Simons et al. 2003 ) or wavelets \Daubechies and 



Pa^/| 19881 Olhede and Walden 200 2j), with which they nevertheless share strong con- 
nections ^lly and Park 1995 ; Shep p and Zhang\2000\\Walter and 5/?g^|2004||2005| ). 



Strict localization of this type remained the driving force behind the development of 
Slepian functions over fixed geographical domains on the surface of the sphere i\Al-\ 
\bertella et <2/.||1999l \Miran ian 2004[ \Simons et al. ||2006] ) — which have numerous 



applications in geodesy ^Ibertella and Sacerdote 



2QQT| \Han eraL\\200Ssi\ \Simons 



I Spatio spectral concentration in the Cartesian plane 



3 



\and Dahlefi\2{){)6), geomagneti sm {Simons et al 2009) Schott and Thehault 20111), 
geophysics \Han and Ditmar\2^{n ^^Han et a l 2008b ; Han and Simons 200 8j ^Harig\ 
\et al.\20lO^, biomed ical { Maniar and Mitrap .005 ; Mitra and Maniar 2006 ) and plan- 
etary {Evans et a/. 1120 10 Han 2008; Han et al. 2009; Wieczorek and Simons 2005) 
science, and cosmology {Dahlen and Simons ^^2QQ'^^ ^Wieczorek and Simons^(lQQl ) 
— as opposed to approaches using spherical wavelets ( Chambodut et al. [2005 Fay 
et al. 2008; Fengler et al. 2007; Freeden and Windheuser 1991 \ Holschneider et al. 
2003 > ^KidoelaL^2W3 ; McEw en et al, 2007 , Panet e t al, 2006, .Schmidt et al.;200^ . 



needlets \Marinucci et al. 2008]), splines ^mirbekya n et a/.||2008[[Za/ et a/.||2009| 
Michel and Wolf 2008 ), radial basis functions ^Freeden and Michel 1999; Schmidt^ 
et a l. 2007 ), coherent states {Hall and Mitchell\.2002, yKowalski and Rembielinski\ 
2000; 7ggm( 2r^|1995||1996| ), or other constructions ^Simons et al.\l997] , which have 
all come of age in these fields also. Finally, we note that the very choice of the cri- 
terion to define localization is open to discussion \Narcowich and Ward\\99^ \Parks 



and Shenoy\l990\\Riedel and Sidorenko\1995\\Saito\2001\\Wei et al.\2010y, in parti c 
ular, it need not necessarily include only quadratic forms {Don oho and Stark\l9S9\ . 
In multiple Cartesian dimensions, in practice: in the plane, the space-frequency 
localization problem has received remarkably little attention beyond the initial treat- 
ment by Slepian himself ^Slepian\\l964\ , who restricted his attention to concentra- 
tion over circular disks in space and spectral space (see also \Brander and DeFacio\ 
|1986t[v^ de Ville'er al. 2002). In this situation, as in the one-dimensional case, the 
concentration operator commutes with a second-order Sturm-Liouville differential 
operator, which greatly facilitates the (numerical) analysis. The scarce recent work 
on two-dimensional Slepian functions has focused on their being amenable to gener- 
alized Gaussian quadrature ^de Villiers et al. |20"Q3} Ma et a l. 1996; Shkolnisky 2007) 
or on using them for multiscale out-of-sample extensions {Coifman and Lafon{2006 ) 
without straying from circular regions of interest. On square or rectangular domains 
Slepian functions formed by outer products of pairs of the pswf \Borcea et al. [2008 



\Hanssen\199 7) have correspondingly square or rectangular concentration regions in 

spectral space, which is undesirable in geophysical applications ( Simons et al. |2003 ). 

Geographical regions are not typically squares, circles or rectangles — although 

usually we will require the spectral support of the Slepian functions to remain disk- 



like to enable isotropic feature extraction \Zhang\\99A ). In those studies where Carte- 
sian domains of arbitrary geometry have been implicit or explicitly considered, as for 
applications in image processing i^Zhou et f2/.jl9 84]), radio -astrono my {^Jackson et al.\ 
|1991| ), medical \Lindquist et al.\2006\ \Walter an d Soleski 2008; Yang et al.\2002\ , 
radar or seismic imaging {Borcea et al. 2008 ) or in the specific context of multi- 



dimensional spectral analysis {Bronez 1988; Liu an d van Veen\\l992^ , the Slepian 
functions are obtained directly from the defining equation, i.e. by numerical diag- 
onalization of the discretized space-bandwidth projection operator. This can be im- 
plemented explicitl y ^Percival and Walden\l993\ or by iterate d filtering and bound- 
ing ^Jackson et a/.||1991 ) as in the Papoulis { Kennedy et a/.||2QQ8t \Papoulis\\l91 5\ 
or Lanczos { Golub and van Loan)\\9'^9 ) algorithms. Given the characteristic step- 
shaped eigenvalue spectrum of the operators {Daubechies\l9^^ 1992[ Slepian\l976 



1983| ) such procedures are not uniformly stable {Bel l et al.\\l993\ , and for general. 



non- symmetric domains, no better-behaved commuting operator exists {Brander and 



E 



Frederik J. Simons, Dong V. Wang| 



DeFacio 1986[ Grilnbaum et al. 1982[ Parle tt and Ww 1984 ). In this case we prefer the 
numerical solution of the Fredholm integral equation {Tricomi\\91Q) with which the 



concentration problem is equivalent, by the Nystrom method { Nystrdm\\9?)Q) , using 
classical Gauss-Legendre integration { Press et al.\\992) . When both the spatial and 
spectral regions of concentration are irregular in shape and of complete generality, 
diagonalization of the projection operator is our only recourse, as we illustrate. 

As we have in a sense come full circle as regards "Slepian's problem", at least 
concerning the scalar case, and because of its deep connections within the framework 
of reproducing-kemel Hilbert spaces ( |Ar6>^5'zq/>z|195Q[|5/m6>^5'|2Q10[|7a6>|1967| ), our 
contribution has the character of a review. We take the reader from the beginnings of 
the theory in one linear dime nsion \Slepian and Pollakfl961) to its general ization on 
the surface of the unit sphere ^Simons andPahlen 2Q06[ Simons etal.\2006), an d back 
to the two Cartesian dimensions of the title of this paper, where, after Slepian\^\96A\ , 
there remained some work to be done. Our discussion is only as technical as required 
to make the theory available for applications in the natural sciences. Focusing on 
the practical and the reproducible, we are distributing all of our Matlab computer 
routines freely over the World Wide Web. 



2 Spatiospectral concentration on the real line 

We use t to denote time or one-dimensional space and uj for angular frequency, and 
adopt a normalization convention {Mall at 199^ ) in which a real- valued time-domain 
signal f{t) and its Fourier transform F{uj) are related by 



/OO />CX) 

-CX) J — OO 



(1) 



Following Slepian, Landau and Pollak ( Landau and Pollak\ 1 96 1 [ Slepian and Pollak 
[1961] ), the strictly bandlimited signal 

nW 

g{t) = {27r)-' / G{uj)e'^'duj, (2) 

J-w 

that is maximally concentrated within \t\ < T is the one that maximizes the ratio 



x = 



\t)dt 



(3) 



g\t)dt 



Bandlimited functions g(i) satisfying problem ^ have spectra G{uj) that satisfy 



w 



D{uj,uj')G{uj')duj' = \G{uj), \uj\ <W, 



w 



D{uj,uj') 



sin[T(co' — uj' 
7t{uj — uo') 



(4a) 
(4b) 
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The corresponding time- or spatial-domain formulation is 



j D{t,t')g{t')dt' = \g{t), tGR, 
sm[W(t - t')] 



D(t,t') = 



(5a) 



(5b) 



TT(t - t') ' 

When the prolate spheroidal wave functions gi{t)^g2{t)^ . . . that solve eq. ^ are 
made orthonormal over |t| < oo they are also orthogonal over the interval \t\ <T\ 



/oo 
-co 

It now follows directly that 

oo 



/ gag(3dt = XaSal3' 



and 



Ealit) 



w 



(6) 



(7) 



A change of variables and a scaling transform eq. Q into the dimensionless eigen- 
value problem 



/■ 



D{x,x')^lj{x')dx' = AV^(x), 



D{x,x') 



^m[TW{x - x')] 



(8a) 



(8b) 



7r(x — x') 

The eigenvalues Ai > A2 > . . . and eigenfunctions ?/^i(x), ?/^2(^), • • • depend only 
upon the time-bandwidth product TW . The sum of the concentration values A relates 
to this as 



N 



ID 



00 ^1 



D{x^x) dx 



2TW _ {2T){2W) 

7T 27T 



(9) 



The spectrum of eq. ([8]) has near-unity and near-zero eigenvalues separated by a nar- 
row transition band (Landau 1965; Slepian and Sonnenblick 1965). Thus, N^^ , the 
"Shannon number", roughly equals the number of significant eigenvalues. In other 
words { Landau and Pollald\l962\ , it is the approximate dimension of the space of 
signals that can be simultaneously well concentrated into finite time and frequency 
intervals. The eigenvalue- weighted sum of the eigenfunctions will be nearly equal to 
the constant in eq. ^ over the region of concentration, as 



N' 



E^«^'w^E^-^'w 



N^^/(2T) if -T <t<T, 
otherwise. 



(10) 



The integral operat or in eq. ^ conmiutes w ith a Sturm-Liouville differential 
operator { Slepian\l9S3\ Slepian and Pollak\l96l\ , to the effect that although x 7^ A, 
we can solve for the functions 7/; also from 



d 
dx 



(l-x^) 



2n# 



dx 



{N 



1D^2^2 



^ = 0, 



\x\ < 1. 



(11) 
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At N discrete values of x = 0, ..., A^ — 1, the eigenfunctions iIj{x) of eq. ( 1 1 ) can be 



found by diagonalization of a simple symmetric tridiagonal matrix ^Grunbaum\l9^l 
Percival and Walden\1993\ \Slepian\l97^ with elements 



[{N- 



2a;)/2]2cos(27rVF), 



'-X x-\-l 



{x^l){N -x-l)/2. 



(12a) 
(12b) 



The matching concentration eigenvalues A can then be obtained directly from eq. ([8]). 
Compared to those, the Sturm-Liouville eigenvalues x are very regularly distributed 



and thus the computation of the Slepian functions via diagonalization of eq. (12) is 
always stable { ^Percival and Walden\1993) . 



3 Spatiospectral concentration on a sphere 

We use r = (^, 0) to denote a location at colatitude and longitude (j) on the unit 
sphere i? = {r : ||r|| = 1}, and adopt a normalization convention (Dahlen and 
Tromp 1998} Edmonds 1996| in which a real- valued time-domain signal /(f) and its 
spherical-harmonic transform fim at degree / and order m are related by 



CX) I » 

/(^) = E E fimyimir)^ fim = / f {r)Yim{r) df2 



(13) 



Following Simons, Wieczorek and Dahlen {Simons et al. 2Q06| ) the strictly band- 
limited signal 

L I 

^(^) = E E 9lmYlm{v) (14) 

that is maximally concentrated within a region i? C i? is the one that maximizes the 
ratio 

f{Y)dQ 



A 



R 



g\v)dQ 



(15) 



Bandlimited functions ^(f ) satisfying problem ( p3] ) have spectra gim that satisfy 

L V 

y^ y^ Dlm,l'm'9l'm' = ^9lm, < I < L, (16a) 



l'=Oi 



Dlmd' 



yimyi'm' dQ. 



(16b) 



Through the addition theorem { Dahlen and Tromp\\99^ , the corresponding spatial- 
domain formulation is 



I D{v,v')g{v')dQ' = \g{v), r e f2, 

Jr 



(17a) 
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D(f,f') = E(^)^'(f-r'), (17b) 

7— n V / 



^=0 



where Pi is the Legendre function. When the functions ^i(r), ^2(r), . . . , g{L-\-iy (r) 
that solve eq. ^Tf) are orthonormal over i? they are also orthogonal over the region R: 



/ 9a9i3 dQ = Saf3, / 9a9(3 dQ = Xa^a^. 



(18) 



From eqs. ([T7])-([T8]) we then immediately obtain the relations 

I)(r,fO= E 9o.{v)g^{v') and ^ ^'(^) = H^' ^^^^ 



Asymptotically, as the spatial area A ^ 0, and L ^ oo, eq. (17) becomes (Si- 
\mons and Dahlen\2001\\Simons et a/.|2006| ) 



/./ 



D«,«')<(«')<i«: = A*(«), (20a) 



27r lie -III 

where the scaled region R^ now has area 47r and Ji is the first-order Bessel function 
of the first kind. As in the one-dimensional case ([9]), the eigenvalues Ai > A2 > . . . 
and eigenfunctions V^i(e), V^2(C)7 • • • depend only upon the product of the maximal 
degree L and the area A. The sum of the concentration values A , the space-bandwidth 
product or "spherical Shannon number", N^^, once again roughly the number of 
significant eigenvalues, is: 



(L+l)^ 

N^''= E ^« = E E ^im,im= I D{r,r)df2, (21) 

q;=1 ^=0 171= — I 



a=l ^=0 m=-l "^^ 

/ I)(|,|)da = (i: + 1)'^. 



The first A^^^ orthogonal eigenfunctions ^a, a = 1,2,..., A^^^, with significant 
eigenvalues Ac^ ^ 1, provide an essentially uniform coverage of the region R, reach- 
ing the constant of eq. ( 19 ) in the following sense: 

E ^'^alir) « E ^"^"(r) - { otherwise. (^2) 

Irrespectively of the particular region of concentration that they were designed for, 
the complete set of bandlimited spatial Slepian eigenfunctions ^i, ^2, • • • , 9{L+iy is a 



basis for bandlimited scalar processes anywhere on the surface of the unit sphere {Si- 



\mons and Dahleri}\2^^6\ \Simons et f2/.||2QQ6] ). The reduced set ^1, ... , 9n^d, with 
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eigenvalues A ^ 1, is an approximate basis for ban dlimited processes that are pri- 
marily localized to the region R {S imons et al.\2009} . Thus N^^ is the approximate 
dimension of the space of signals that can be simultaneously well concentrated into 



finite spatial and spectral intervals on the surface of the unit sphere. Eq. ( [2Q| ) depends 
only on this combination of bandwidth and spatial area, as in the one-dimensional 
case, eq. ([5]). 

When the region of concentration is a circularly symmetric cap of colatitudinal 
radius O, centered on the North Pole, the colatitudinal parts g{6) of the separable 
solutions to eq. ([17]), 



^(^,</>) = < 



' V2g(l9)cos(m^) 

m 

V^g(l9)sin(m^) 



if -L < m < 0, 

if m = 0, 

if < m < L, 



(23) 



are, owing to a commutation relation \GrUnbaum et al.\l9^2) , identical to those of 
the Sturm-Liouville equation 



d 
dji 



[jl — COS0)(1 — jl ) 



2^^g 



dfi 



(24) 



X + L(L + 2)/i- 



m^{li — COS0) 



g = 0, 



with /i = cos and x 7^ A. At constant m, their values gim in the expansion ( 14 ) can 



be found by diagonalization of a simple symmetric tridiagonal matrix ( Grunbaum 
\et al.\l9'S2\ \Simons et al. |2006| ) with elements 



T,, = -;(; + !) cos 0, 




(25a) 


T„+i= [l{l + 2)-L{L + 2)\^ 


1 {I + 1)2 - m^ 
{2l + l){2l + 3)' 


(25b) 



When the region of concentration is a pair of axisymmetric polar caps of common 
radius O centered on the North and South Pole, the g(^) solve the Sturm-Liouville 
equation 



X + Lp{Lp 



(26) 



3)^2 



m^{li^ — cos^ O) 



= 0, 



where Lp = L or Lp 



L — 1 depending whether the order m of the functions giO^cj)) 
in eq. j23] ) is odd or even, and whether the bandwidth L is odd or even ( Grunbaum\ 
et a/.|[l98 2i Simo ns and Dahlen\2006) . The expansion coefficients of the optimally 
concentrated antipodal polar-cap eigenfunctions require the numerical diagonaliza- 
tion of the symmetric tridiagonal matrix ( Simons and Dahlen\2{){)6 ) with elements 
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Ti; = -i{i^i)cos^e- 



2/ + 3 
[(/-2)(/ + l)-L^(L^ + 3)] 



[(/ + 1)2 - m^] (27a) 

1 2 Sm^ - /(/ + 1) 



rpp 



21 



3 3 (2/ + 3)(2/-l) 



[/(/ + 3)-L^(L^ + 3)] /[(/ + 2) 



m 



{i + if 



m 



(2/ + 5)(2/ + l) 



. (27b) 



Every other degree in the expansion for the equatorially (anti-)synimetric double- 
polar cap functions is skipped \Simons and Dahlen\2006 ), hence the / + 2 subscript 
for the elements off the main diagonal in eq. (27b ). The concentration values A can be 
determined from the defining equations (p^ or (|17]). Compared to these, the Sturm- 
Liouville eigenvalues x in eqs. ( [24| ) and ( [26| ) are very regularly spaced, and thus the 
computations for these special cases are inherently stable. 



4 Spatiospectral concentration in the Cartesian plane 



We now turn to the multidimensional Cartesian case, first discussed by |5/g/7/a^ ( |1964| ), 
noting that we have set ourselves up for a result that is analogous to eq. pO|). Indeed, 
in the asymptotic regime of an infinitely large bandwidth and an infinitesimally small 
region on the surface of the unit sphere, the spherical and Cartesian concentration 
problems are of course equivalent ( Simons and Dahlen\2001\\Simons et aL\2006 ). 

Another note concerns the equivalence between the temporal or spatial and spec- 
tral forms of the concentration problems. In writing eqs. Q-ijS]), and eqs. (p^-(p7|), 
respectively, we have exclusively considered strictly bandlimited, time- or space- 
concentrated Slepian functions. Strictly time- or spacelimited, band-concentrated func 
tions can be obtained via an appropriate restriction of the integration domains in these 
equations, and the resulting new functions can be obtained from the old ones by sim- 
ple truncation and rescaling. Both for the one-dimensional and spherical situations 
this distinction is usually made more explicitly elsewhere ^Landau and Pollak\l96l 



Simons et al. ||2Q06{ \Slepian a nd Pollaid\\96l) , and in what follows we once again 
treat both cases separately. 

For applications in the geosciences, we focus on the two-dimensional, flat geom- 
etry of geographical maps: the Cartesian plane R^, with a spatial coordinate vector x 
and a spectral coordinate vector k. To wit, we formulate the concentration problem as 
follows: in spite of the Paley -Wiener theorem, which states that functions cannot be 
spatially and spectrally restricted at the same time \Daubechies\ 1 992[ Ma/to|1998] ), 
can we construct (real) functions that are localized to, say, the shape of Belgium, in 
x-space, while having a Fourier transform localized to, say, the (Hermitian) shape of 
a pair of triangles, in k- space? Yes we can! 
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4.1 Preliminary considerations 

A real- valued, square-integrable function /(x) defined in the plane has the two- 
dimensional Fourier representation 



/(x) = (27r)-2 / F(k)e^^-^dk, F{k) = f /(x)e-^^-^6ix, 

with F(k) = F*(— k). The Fourier orthonormality relation is 

(2^)-' / e'^-("-"')dk = (5(x,xO, 



(28) 



(29) 



which defines the delta function in the usual distributional sense 

/ /(x0(^(x,x06^x' = /(x). (30) 

Likewise, in the spectral domain we may write 

(2^)-' / e'(^-^'>"^x = (5(k,kO. (31) 

By Parseval's relation the energies in the spatial and spectral domains are identical: 

/ /2(x)dx=(2^)-2/ \F{k)\'dk. (32) 

JR2 J^2 

4.2 Spatially concentrated bandlimited functions 

We use ^(x) to denote a real function that is bandlimited to /C, an arbitrary subregion 
of spectral space, 

^(x) = (2^)-2/G(k)e^^-^k, (33) 



Following Slepian ( 1964| ), we seek to concentrate the power of ^(x) into a finite 



spatial region IZ C [R^ of area A, by maximizing the energy ratio 

/ ^^(x)(ix 
Jn 

/ ^^(x)(ix 

JR2 



(34) 



Upon inserting the representation (33 ) of ^(x) into eq. (34) we can express the con 



centration in the form of the Rayleigh quotient 

/ f G%k)D{k,k:)G{k:)dkdk: 



(35) 



\G{k)\'dk 
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where we have used Parseval's relation ( 32 ) and defined the positive-definite quantity 

D{k, kO = (27r)-2 / e^(^'-^>^ ^x, (36) 

Jn 

which is Hermitian, D(k, k') = D*(k',k). BandHmited functions ^(x) that maxi- 



mize eqs. (|34|)-(|35]) solve the Fourier-domain Fredholm integral equation 

L>(k, kO G{k') dk' = AG(k), k G /C. (37) 



Jk 



Comparison of eq. ( |36| ) with eq. p\} leads to the interpretation of the spectral-domain 
kernel D{k, k') as a spacelimited spectral delta function. We rank order the concen- 
tration eigenvalues so that l>Ai>A2>...>0. Upon multiplication of eq. ( [37] ) 
with e*^"^ and integrating over all k G /C, we obtain the corresponding problem in 
the spatial domain as 



/ D(x, xO ^(xO ^x' = A^(x), X G [R^ 
Jn 

L)(x,xO = (27r)-2 / e^^-(^-^') dk. 



(38a) 



(38b) 



Comparison of eq. ( |38b| ) with eq. ([29]) shows that the Hermitian spatial-domain ker- 
nel I)(x, x') = I)*(x',x) is a bandlimited spatial delta function. The bandlimited 
spatial-domain eigenfunctions gi (x) , ^2 (x) , . . . may be chosen to be orthonormal 
over the whole plane R^, in which case they are also orthogonal over the region 7Z: 



/ gag (3 d-X. = Sa(3, / QaOfS <^X = XaSa(3- 

JR2 J-j^ 



(39) 



In this normalization the eigenfunctions of eq. ( 38 ) represent the kernel in ( 38b ) as 

CX) 

Z)(x,x') = ^5a(x)5a(x')- (40) 



This form of Mercer's theorem \Flandrin\ 1 998[ Tricomi\ 1 970 ) is verified by substi- 
tuting the right hand side of eq. (|40|) into eq. ( |38al ) and using the orthogonality ( [39] ). 



4.3 Spectrally concentrated spacelimited functions 

We use /i(x) to denote a function that is spacelimited, i.e. vanishes outside the arbi- 
trary region IZ of physical space: 



ff(k) 



h{x.)e' 



ikx 



dx. 



(41) 



n 
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To concentrate the energy of h{x.) into the finite spectral region /C, we maximize 



\H{k)\^dk 
\H{k)\'^dk 



(42) 



Upon using eq. (32) and inserting the representation (41 ) of ilf (k) into eq. (42) we 
can rewrite A in the form 



JnJn 



Jn 



dx 



(43) 



where we again encounter the quantity (38b), 



I)(x,xO = (2^)-' / e'^-^"-"') ^k. 

JK 



(44) 



Once again by Rayleigh's principle, spacelimited functions /i(x) that maximize the 



quotient A in eqs. (42)-(43 ) solve the spatial-domain Fredholm integral equation 



L' 



Di(K,^) h{^) d^ = A/i(x), X G 7^. 



(45) 



Eq. ([45]) is identical to eq. (38 ) save for the restriction to the domain 1Z. The eigen- 



functions /i(x) that maximize the spectral norm ratio (42) are identical, within the 



region 1Z, to the eigenfunctions ^(x) that maximize the spatial norm ratio (34). The 



associated eigenvalues 1 > Ai > A2 > . . . > are a measure both of the spatial 
concentration of ^(x) within the region 7Z and of the spectral concentration of /i(x) 
to the wave vectors k G /C. Identifying 



/i(x) 



_/<7(x)ifxe7e, 
otherwise, 



(46) 



the normalization is such that 



R2 Jn 



(47) 



The null- space consisting of all non-bandlimited functions that have both no energy 
outside IZ nor inside /C is of little consequence to us. 
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5 Slepian Symmetry 



Under what has been called "Slepian symmetry" {Brander and DeFacio\\9'^6) the 



solutions to eqs. ( [38] ) and ( [45] ) take on particularly attractive analytic forms. Antic- 
ipating a switch to polar coordinates x = {r^6) we introduce Jm{k), the Bessel 
function of the first kind and of integer order m {\bramowitz and Stegun 1965 



Gradshteyn and Ryzhik\2006 ). The Bessel functions satisfy the symmetry condition 
J-m{k) = (-l)^J^(A:), the relation 

oo 

l = j2(fc) + 2^ J^(fc) (48) 

and the identity ^Jeffreys and Jeffrey s\l9^^ ) 

CX) 

Jo(fc||x - x'll) = Jo(fcr) Jo(fcr') + 2 ^ Jm{kr)Jm{kr') cos[m{e - 0')]. (49) 

171=1 

Furthermore, we have the particular formulas 

r27r 



Jo{k) = (2^)-^ 



Ak c 







'd(9, Ji{k) = k-^ jQ{k')k'dk', (50) 
io 



Ji/2{k) 



sin/c, 



the derivative identity -^[kJi{k)] = kJo{k) and the limits 



fc^o k 



and 



Jm\k) 

lim — — > 0. 



(51) 



(52) 



5.1 Circular bandlimitation 

Limitation to the disk-shaped /C = {k : ||k|| < K} allows us to rewrite the spatial 
kernel of eq. (44) using polar coordinates as follows. Letting k = ||k|| = (k • k)^/^ 
and 6 be the angle between the wave vector k and the vector x — x^ 



k-(x-xO = fc||x- x'll cos i9, 



and using eq. (50), we obtain an expression alternative to eq. (44) as 

rK /.27T 



I)(x,xO = (2^)"' / / "e'^ll"-"'ll 
JO io 



dOkdk, 



K 



(27r)-^ / Jo(/c||x-x'||)/cdfe, 



KJi(K||x-x'| 
2^||x-x'|| 



(53) 

(54a) 
(54b) 
(54c) 
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We notice the symmetry I^(x, x') = D{x.'^x.) = I)(||x — x'||) and the equivalence 
of eq. ( 54c ) with the asymptotic expression (|20b[), as expected ^Simons and Dahlen 



2007} \Simons et f2/.||2QQ6"] ). With \Coifman and Lafori\ ( |20Q6| ) and through eq. (|5T])^ 



we furthermore note that both eqs. ( |54c| ) and (5b) are n-dimensional versions of the 
general form [i^/(27r)]^/2 Jn/2(i^||x - x'||)/(||x - x'H^/^), for n = 2, and n = 1 
(where K = W), respectively. As eq. ([52]) now shows, 



I)(x,x) = ^, (55) 



and thus eqs. ( 39 )-(40 ) allow us to write the sum of the concentration values A as the 



"planar Shannon number" 



Ar2D = y^ A« = / D(x,x) rfx = K^ 

1 JTl 



47r (27r)2 



(56) 



where A is the area of the spatial region of concentration TZ. As in the one-dimensional 
case, eq. ([9]), the Shannon number is equal to the area of the function in the spatiospec- 
tral phase plane times the "Nyquist density", as expected { Daubechies\l9SS\ 1992| ). 



Roughly equal to the number of significant eigenvalues, N^^ again is the effective 
dimension of the space of "essentially" space- and bandlimited functions in which 
the reduced set of two-dimensional functions ^i , ^2 , • • • , 9n'^^ ^^Y act as an efficient 



orthogonal basis. Using eqs. ( [55] ) with ( [4Q| ) we furthermore see that the sum of the 
squares of all of the bandlimited eigenfunctions, independent of position x in the 
plane, is the constant 

q;=1 

Likewise, since the first N'^^ eigenfunctions ^1 , ^2 , • • • , 9n^^ have eigenvalues near 
unity and lie mostly within IZ, and the remainder gN'^D-\-i^9N'^'^-\-2^ • • • ^9 00 have 
eigenvalues near zero and lie mostly outside IZ in R^, we expect the eigenvalue- 
weighted sum of squares to be 



00 iV^D 

Yl ^^9l{^) ^ Yl ^oc9l{^) ^ \ 



/V-2 

'N^^/Aif^€Tl, 
otherwise. 



This heuristic finding is of great importance in the analysis and representation of 
signals using the Slepian functions as an approximate basis, as much as for their use 
as tapers to perform spectral analysis on data from which we thus expect to extract 
all relevant statistical information with minimal loss or leakage near the edges of the 
region under consideration ( |Tya/Jg^|1990| ). 



5.2 Scaling analysis I 

Introducing scaled independent and dependent variables 



Spatiospectral concentration in the Cartesian plane 



15 



we can rewrite eqs. ( |38) ) and ( |54c| ) as 



L} 



with 7^*, of area 47r, the image of the concentration region IZ under ( 59 ), and 



D(.t^) 






(60a) 



(60b) 



Eq. ( [60| ) reveals that the eigenvalues Ai, A2, . . . and the scaled eigenfunctions V^i(C), 
^^2(^)5 • • • depend on the maximum circular bandwidth K and the spatial concentra- 
tion area A only through the planar Shannon number N^^ , as they do in the one- 



dimensional, eq. ([8]), and asymptotic spherical, eq. ( [20| ), cases. Under this scaling the 
Shannon number remains 

^2D 



N' 



2D 






(61) 



5.3 Circular spacelimitation 

If in addition to the circular spectral limitation, physical space is also circularly lim- 
ited, in other words, if the spatial region of concentration or limitation 7^ is a circle 
of radius R, then a polar coordinate, x = (r, 0), representation 



9{r,e)={ 



>/2g(r) cos(m6>) if m < 0, 
g(r) if m = 0, 

>/2g(r) sin(m6>) if m > 0, 



(62) 



may be used to decompose eq. ( [6Q| ) into an infinite series of non-degenerate fixed- 
order eigenvalue problems, one for each order ±.m. We do this most easily by insert- 
ing eq. ( [49] ) into eq. ( |54b| ), the obtained result with eq. ( [62] ) into eq. ( |38a| ), and using 
the orthogonality of the functions . . . , >/2 cos(m^), . . . , 1, . . . , >/2 sin(m6>), . . . over 
the polar angles < 6> < 27r. The integral equation for the fixed-m radial func- 
tion g(r) is then given by 

D(r, r) g(r') rdr = Ag(r), (63a) 



/ 

Jo 



with the fixed-m symmetric kernel 

D{r,r')=K^j Jm{Kpr) Jm{Kpr')pdp. 
Jo 



(63b) 



We rank order the distinct but pairwise occurring eigenvalues obtained by solving 
each of the fixed-order eigenvalue problems ([63)) so that l>Ai>A2>--->0, 



and we orthonormalize the associated eigenvectors gi, g2, . . . as in eq. (39 ) so that 

27r / ga{r)g(3{r) rdr = 6^(3, 27r / gc(r)g/3(r) rdr = XaSafS- (64) 



We shall denote the radial part of the functions /i(r, 0) by h(r) in Sect. 5.6 
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5.4 Scaling analysis II 

Finally, the scaling transformations 

i = r/R, i' = r'/R, ^(0 = g(r), ^(O = g(^'), 

convert eq. ([63]) into the scaled eigenvalue problem 



Jo 



AV^(0, 



(65) 



(66a) 



with the fixed-m kernel 



D{^, O = 47V2^ / Jm {2Vn^p^) Jm {2VN^pe) pdp, (66b) 

Jo 



which is dependent only upon the Shannon number 



N' 



2D 



K' 



(67) 



The number of significant eigenvalues per angular order m is ^Simons and Dahlen 
\2mTWmons et al\2m6) 



N. 



2D 



= / D(C,0?rf? = 4iV2i^ / / J^(2ViV^K)p*CC (68) 
Jo Jo Jo 



2N 



~{2m + l)VN^J„,{2VN^)Jm+i{2VN^ 

m 

1 - Jo'(2Va^) - 2 ^ J^(2VAr 



(69) 



m 



r2D ^ 



n=l 



The complete Shannon number N'^^ is preserved, inasmuch as, per eqs. ( 68 ) and (48 ), 

^2D ^ ^2D ^ 2 V iV^^ = 4iV2D / / pdp^di = iV2D. (70) 

Jo Jo 



m=l 



5.5 Sturm-Liouville character and tridiagonal matrix formulation 



As noted by |5/g/7/a^ ( |1964| ), eq. ([66]) is an iterated version of the equivalent "square- 
root" equation 



2Vn^ [ jm{2VN^^e)^{e)ede = ^^{0' 

Jo 



(71) 



To see this, it suffices to substitute for V^(^') on the left hand side of eq. (71 ) using 



eq. (71 ) itself. Thus, the eigenvalues Ai, A2, . . . and eigenfunctions ^l^i{^)^ il;2{C)i- • • 
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of eq. ( 66 ) may alternatively be found by solving the equivalent equation ( 71 ). A fur- 
ther reduction can be obtained by substituting a scaled Shannon number or bandwidth 
and rescaling, 



to yield the more symmetric form 



(72) 



(73) 



On the domain < ^ < 1 the (/:?(^) also solve a Sturm-Liouville equation, 



d 



(1-r 



dip 



1/4 - m^ 



e 



m'^'e] ^ = 0, lei < 1, (74) 

) reduces to the one-dimensional equa- 
tion ( jj_), as can be seen by comparing eqs. (67]) and ^ after making the identifica- 
tions K = W and R = T. The fixed-m Slepian functions (^(^) can be determined by 
writing them as the infinite series 



for some x 7^ A. When m = ±1/2 eq. d74 



^(0 = m!r+^E 



dl^' omO/i ot2 



^=0 



(/ + m)! 



Pr(l-2^^), 0<^<1, 



(75) 



whereby the P[^^{x) are the Jacobi polynomials ^bramowitz and Stegun 1965 ): 



;^(n + m)!(^-n)! 



!2n!' 



(76) 



By extension to ^ > 1 they can also be determined from the rapidly converging 
Bessel series 



^{0 



E 



m! 



7 ^ (/ + m)! 



G^^/! Jm+2i+i{ci) 



ci 



, ^gR+. 



(77) 



From eq. (52) we confirm that (p{0) = 0. In both cases, the required fixed-m ex 
pansion coefficients di can be determined by recursion { Bouwkamp^\9Al Slepian 



1964| ), which is, however, rarely stable. It is instead more practical to determine them 



as the eigenvectors of the non-symmetric tridiagonal matrix {de Villiers et al 2003 
that: 



\Shkolnisky ^2007) that is the spectral form, with the same eigenvalues, of eq. ( [74] ), 

c^(m + / + l)^ 



(2/ + m + l)(2/ + m + 2)' 



Til 



2/ + m- 



1 



2/ + m- 



c 
~2 



m 



t-ii+i 



(2/ + m)(2/ + m + 2) 
c' {I + If 



(2/ + m + 2)(2/ + m + 3)' 



(78a) 

(78b) 
(78c) 
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where the parameter / ranges from to some large value that ensures convergence. 
Finally, the desired concentration eigenvalues A can subsequently be obtained by 
direct integration of eq. ( [60| , or, alternatively ^Slepian\1964\ , from 



A = 2-f^VN^, with 



^+1/2^0 



2^+i(m + l)! 



Kl=0 / 



(79) 



Neither procedure is particularly accurate for small values A, though these will be 
rarely needed in applications. A more uniformly valid stable recursive scheme is de- 
scribed elsewhere ( Slepian and Pollak\ 1 96 1 [ Xiao et al.\2001\ . The overall accuracy 
of the procedures chosen can be verified with the aid of the exact formula eq. (69), as 
for the fixed-m eigenvalues, A^^^ = X^^i A^. 



2- 

II 0- 

E 

-2- 


A, = 1 .000000 1 " 




— 1 



-2 



U= 0.999998 



o-A- 



U= 0.999974 



_/V. 



II 

E 



R= 0.999738 



./V 



V 



Ar 



a = 

1 


2 


. \n 


0.999997 1 - 




v 


1 



a = 3 



V 



'\l^ 0.998757] 



k 



Pi= 0.898354 



U= 0.999930 



Vv- 



U= 0.986938 



A/\/^^^^^ 



Pi= 0.609977 



U= 0.999121 



-/\/U- 



X= 0.915710 



-^\/\:/\y\^r^ 



X= 0.245783 



U= 0.992470 



Ar 



-At- 



\k^ 0.955293 




01 501 501 501 5 

radius ^=r/R radius ^=r/R radius 5=r/R radius ^=r/R 



Fig. 1 Radial dependence of the first four eigenfunctions of eq. J63}, the bandlimited Slepian functions 
ga (r), a = 1,2,3,4, of fixed order m = (top) to m = 4 (bottom), calculated from eq. f/V) truncated at 
/ = 84. The Shannon number A^^^ = 42 and the bandwidth c ~ 13. Black curves show the concentration 
within the scaled spatial interval < ^ < 1 and grey curves show the leakage into the rest of the real line, 
which is truncated at ^ = 5. Labels show the eigenvalues, Xa, calculated from eq. f79) . These express the 
quality of the spatial concentration. Spacelimited, spectral-domain Slepian functions are shown in Fig. [2] 
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5.6 Numerical examples 



With the eigenfunctions (p{^) calculated, by series expansio n ^de Villiers et f2/ .|2QQ3t 
\Shkolnisky\20Qn) via eqs. ( [75] ) o r ([77] ) or even by quadrature \Zhang\ 1 994] ) of eq. ( [73] ), 
we rescale them following eqs. (|72[) and (65 ) to the Slepian functions ga(r) normal- 



ized according to ( [64] ). The sign is arbitrary: the concentration criterion is a quadratic. 

The four most optimally concentrated eigenfunctions gi(r), g2(r), g3(r), g4(r), 
for the orders < m < 4 are plotted in Fig.[T| scaled for convenience. The associated 
eigenvalues Ai, A2, A3, A4 are listed to six-figure accuracy. The planar Shannon num- 
ber in this example is N'^^ = 42, hence c ^ 13. It is useful to compare the behavior 
of these solutions to those of the spherical concentration pr oblem ^ mons et ah '2006 
Fig. 5.1), with which they are asymptotically self-similar ^Simons and Dahlen2001 
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Fig. 2 Squared Fourier coefficients |if (/c)p of the first four scaled spacelimited eigenfunctions hQ,(r), 
OL = 1, 2, 3, 4, of fixed order m = (top) to m = 4 (bottom), calculated by discrete Fourier transforma- 
tion of the functions shown by the thick black lines in Fig. IT] The Shannon number is A^^^ = 42 and the 
scaled bandwidth c ~ 13. Black curves show the power within the scaled wavenumber interval < k < 1 
and grey curves show the power leaked to the rest of the wavenumber axis, which is truncated at ^c = 5. 
Values of \H(k)\'^ are in decibel (dB), zero at the individual maxima. Labels show the eigenvalues Xa, 
which express the quality of the spectral concentration. The corresponding bandlimited, spatial-domain 
eigenfunctions are shown in Fig.Q] 
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Simons etal 2006). The first zeroth-order "zonal" (m = 0) eigenfunction, gi(r), has 
no nodes within the "cap" of radius R', the second, g2(^), has one node, and so on. 
The non-zonal (m > 0) eigenfunctions all vanish at the origin. The first three m = 0, 
m = 1 and m = 2, and the first two m = 3 and m = 4 eigenfunctions are very 
well concentrated (A > 0.9). The fourth m = 3 and m = 4 eigenfunctions exhibit 
significant leakage (A < 0.1). 

The squared Fourier coefficients \H{k)\'^ of the four best concentrated space- 
limited eigenfunctions hi(r), h2(r), h3(r), h4(r) for < m < 4 are plotted ver- 
sus scaled wavenumber k = k/c, on a decibel scale, in Fig. [2| Shannon number, 
bandwidth, and layout are the same as in Fig. [T] See QlsQwhQro JSimons et aL\2006 
Fig. 5.2) to compare with the spherical case. 

Once a number of sequences of fixed-order eigenvalues has been found, they 
can be resorted to an overall mixed-order ranking. The total number of significant 
eigenvalues is then given by eq. (jTOb. In Fig. |3] we show the mixed-m eigenvalue 
spectra for the Shannon numbers N^ = 3, 11, 24 and 42. These show the familiar 
step shape ( Landau\l965\ Slepian and Sonnenblick\l965 ), with significant (A ^ 1) 
and insignificant (A ^ 0) eigenvalues separated by a narrow transition band. The 
Shannon numbers roughly separate the reasonably well concentrated eigensolutions 
(A > 0.5) from the more poorly concentrated ones (A < 0.5) in all four cases. To 
compare with the spherical case, see elsewhere ( Simons et a;/. [2006 Fig. 5.3). 
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rank 
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Fig. 3 Eigenvalue spectra {\ol versus rank a) of the Cartesian concentration problem, at various Shan- 
non numbers A^^^, for disk-shaped regions, calculated via eq. {79) . The total number of eigenvalues 
is unbounded; only Ai through Aeo are shown. Different symbols are used to plot the various orders 
— 11 < m < 11; juxtaposed identical symbols are ±?7z doublets. Vertical grid lines and top labels specify 
the rounded Shannon numbers A^^^ = 3, 11, 24, 42. 
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Fig. 4 Bandlimited eigenfunctions g(r, 0) that are optimally concentrated within a disk of radius R = 1. 
Dashed circle denotes the region boundary. The Shannon number A/"^^ = 42 and the scaled bandwidth 
c ^ 13. The eigenvalues Xa have been sorted into a mixed-order ranking, o; = 1 — )► 30, with the best 
concentrated eigenfunction plotted on the top left and the 30th best on the lower right. Blue (or white in 
the print version) is positive and red (black in print) is negative; regions in which the absolute value is less 
than one hundredth of the maximum value on the domain are left white (grey at 50% in the print version). 



Finally, Fig. [5] shows a polar plot of the first 30 eigenfunctions ^(r, <p) concen- 
trated within a radius R = 1, defined by equations ( [62| ) using the g(r) shown in 
Fig.[l] The Shannon number is N'^^ = 42, as in Fig.jipThe eigenvalue ranking is 
mixed-order, as in Fig.[3| and all degenerate >/2cos(m^), ^/2sm{m(j)) doublets are 
shown. The concentration factors 1 > A > 0.8983 and orders m of each eigenfunc- 
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tion are indicated. Blue and red colors (in the online version; in print these have been 
converted to grey scales between white and black) represent positive and negative val- 
ues, respectively; however, all signs could be reversed without violating the quadratic 
concentration criteria ( [34] ) and ( [42] ). Elsewhere ( Simons et a/.||20Q6[ Fig. 5.4) these 
can be compared with the spherical case. Other calculation methods (see Sect.[6]be- 
low) may yield slight differences between what should be pairs of eigenvalues for 
each non-azimuthally symmetric and thus fixed non-zero order eigenf unction. The 
mismatch is confined to the last two quoted digits in a typical calculation (compare 
5/m6>^^|2Q10| Fig. 2). 



6 Computational Considerations on General Domains 

While the framework of Sect. |4]is completely nonrestrictive as to the geometry of the 
spectral or spatial concentration regions /C or IZ, the literature treatment has remained 
focused on the cases with special symmetry that we discussed in Sect. |5j for which 
analytic solutions could be found. In what follows we lift all such restrictions and, 
numerically, solve the problem of finding the unknown planar function /(x) that is 
bandlimited/band-concentrated to /C and space-concentrated/spacelimited to 7^, by 
satisfying the Fredholm eigenvalue equation 



/.■ 



I^(x,xO/(xOc^x^ = A/(x), (80a) 



L)(x,xO = (27r)-2 / e^^-(^-^') dk, 



(80b) 
over the domains x G R^ (in which case we retrieve the bandlimited Slepian func- 



tions g of Sect. 4.2 ) or X G 7^ (when we obtain the spacelimited h of Sect. 4.3 ). 



6.1 Isotropic spectral response 



Despite the appealing generality of eq. (80), applications in the geophysical or plan 



etary sciences, e.g. in problems of spectral analysis of data collected on bounded 
geographical domains, is well served by Slepian functions whose spectral support is 
circularly isotropic. As we saw in Sect. |5.1| the integral kernel is then 

_ i^Ji(i^llx-x-ll) 
^ '^" 27r||x-x'|| ■ ^^'^ 

The standard solution to solve homogeneous Fredholm integral equations of the sec- 
ond kind is by the Nystrom method { Nystrdm\l9?)Q\ Press et al.\l99 2\. 



In one dimension, let us write eq. dSOj) in a form reminiscent of the special case 
that we encountered previously in eq. ( |63| ), to which it also applies explicitly, namely 

h 

D{x, x') f{x') dx' = Xf{x). (82) 
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The left-hand side of this equation is to be discretized by a quadrature rule; this 
involves choosing J weights Wj and abscissas Xj such that, to the desired accuracy, 



/ D{x,x')f{x')dx' = ^WjD 



(X, Xj) j[Xj)^ 



(83) 



using which we rewrite eq. ( |82| ) as 

J 



^WjD{xi,Xj)f{xj) = \f{xi), 



(84) 



this time evaluating the unknown right hand side f{x) at the quadrature points as 
well. Written in matrix form, we may solve 



D W f = Af, 



(85) 



identifying the elements of the square kernel Dij = D{xi^ Xj), the samples of the 
unknown function as fj = f{xj), and the entries of the diagonal weight matrix 



Wij = Wj5ij. 
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Fig. 5 Bandlimited eigenfunctions gi,g2, . . . ,95, calculated by eq. J88} using Gauss-Legendre integra- 
tion with 32 node points in both directions, optimally concentrated within the Colorado Plateaus, centered 
on 109.95°W 37.01°N (near Mexican Hat, Utah) of area A ^ 334 x 10^ km^. The concentration factors 
Ai , A2 , . . . , A5 are indicated; the Shannon number is A^^^ = 10. The top row renders the eigenfunctions 
in space on a grid with 5 km resolution in both directions, with the reversible convention that positive 
values are blue (white in print) and negative values red (black in print). The spatial concentration region is 
outlined in black. The bottom row shows the squared Fourier coefficients |G(k)p as calculated from the 
functions gf(x) shown, on a wavenumber scale that is expressed as a fraction of the Nyquist wavenumber. 
The spectral limitation region is shown by the black circle at wavenumber K = 0.0194 rad/km. All areas 
for which the absolute value of the functions plotted is less than one hundredth of the maximum value 
attained over the domain are left white (grey at 50% in the print version). 
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Restoring the symmetry of the problem contained in the symmetry D = D^, and 
provided the weights Wj are positive, we had rather solve 



W D W f = Af, 



(86) 



using Wij = yjwibi^ and f = W • f or indeed fi = ^/wlfi. Having solved this 
system for the unknown function /(x^) at the integration nodes, we can subsequently 
produce either g{x) (anywhere in one-dimensional space) or h{x) (inside the region) 
using eqs. ([82|)-([83]), to the same accuracy. 

In two dimensions, the vectors x = (x, y) and x' = {x\ y') denote points occu- 



pying the plane and eq. ( 80 ) is rewritten as 

K L 



/ i:>(x, x') /(x') ^x' = ^ ^ WkWiD{yi] Xk, yi) f{xk, yi), 

'^^ h — -{ 7 — 1 



(87) 



fc=l ^=1 



where it is implied that wi = wi (x^) and yi = yi {x^). By wrapping the indices k and 
/ on the scalar elements x^ and y^ into the vectorial indices i and j used for x and x' 
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Fig. 6 Cumulative eigenvalue-weighted energy of the first A^^^/2, A^^^ and 2N^^ eigenfunctions that 
are optimally concentrated within the Colorado Plateaus, with the Shannon number N"^^ = 10. The top 
row is in the spatial domain with the concentration region outlined in black. Darkest blue (white in print) 
on color bar corresponds to the expected value {58) of the sum, as shown. The bottom row is in the spectral 
domain with a wavenumber scale that is expressed as a fraction of the Nyquist wavenumber. The spectral 
limitation region is shown by the black circle. Areas with absolute value less than one hundredth of the 
maximum value attained over the domain are left white (grey at 50% in the print version). 
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we arrive at 

J 

^^p(x,;x,)/(x,) = A/(x,). (88) 



This is identical to eq. ( 86 ) as long as we remember that the weights w'^ are pairwise 
products of the one-dimensional WkWi and identifying Dij = D{^i] Xj). 

As to the choice of quadrature, a finely meshed Riemann sum might suffice for 



some applications {Saito 2007 Zhang |1994|), though we shall prefer the classical 



p994]j, 



Gauss-Legendre algorithm {P ress et al.\l992) for increased accuracy at faster com- 
putation speeds. Other options are available as well \Ramesh and Le an\l99\) . The 



analytical formulas for the special cases, e.g. eqs. ( [75] ), (77 ) and (79 ) can be used for 



verification purposes (see the discussion of Fig.|4]in the text), as we have 

Figures [5] and [6] were computed using eq. ( [88] ) with Gauss-Legendre quadrature, 
starting from a splined boundary of the Colorado Plateaus, a physiographic region 
in the United States. Fig. [5] shows the first five bandlimited Slepian functions, gi -^ 
^5, with their eigenvalues, Ai -^ A5. The space-domain functions ^(x), shown in 
the top row, are increasingly oscillatory as their concentration values decrease. At 
the same time their periodograms |G(k) p, shown in the bottom row, remain strictly 
confined within the chosen isotropic spectral region of bandwidth K. All together, 
the Slepian functions uniformly occupy their respective domains of concentration or 
limitation. Fig. [6] shows the eigenvalue- weighted partial sums of the energy of the 
eigenfunctions, in space (top row) and spectral space (bottom row). The progressive 
covering behavior, consistent with eq. ( [58] ), lies at the basis of the success of the 
Slepian functions in being used as data tapers for spectral analysis (i Bronez\l9^mLiu 



and van Veen 1992 Percival and Walden 1993): orthogonal over the domain that they 



cover, smoothly but rapidly decaying to zero at the boundary of interest, and with a 
finite and isotropic spectral response. 



6.2 Arbitrary spectral response 

We finally treat the situation in which both the spatial and spectral concentration 
domains are of arbitrary description — but with Hermitian spectral symmetry, see 
eqs. ([28|)-([29|), if the Slepian functions are to be real- valued. No analytical results can 
be expected for such a case, and we benefit from writing the concentration problem in 
its most abstract form. Writing Q for the operator that acts on a spatial function /(x) 
to return its Fourier transform F(k), we introduce the spatial projection operator 

and the spectral projection operator 

CF{k) = I f ^) '[^ ^ ^' (90) 

^ ^ [0 otherwise. 
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We rewrite the variational equations ( 34 ) and ( 42 ) in inner-product notation as 



A = ^ -. -. -^ = / ^^ . ^^ .V = maximum. (91) 

{Q-^CF.Q-^CF) (QVf.QVf) 

The associated spectral-domain and spatial-domain eigenvalue equations are 

jCQVQ-^jC{CF) = X{jCF), VQ-^CQViVf) = X{Vf), (92) 

where we have made use of the fact that V^ = V and C^ = C are self- adjoint, 



and that Q and Q ^ are each other's adjoints, provided we rewrite eq. (28 ) as a uni 



tary transform. In our notation, the solutions CF yield the Fourier transforms of the 



bandlimited functions ^(x) of Sect. 4.2 while the Vf represent the spacelimited 
functions h{x.) of Sect. 4.3 

For computation in the discrete-discrete case, V and C will simply be matrices of 
ones and zeros, and Q and Q~^ will be the discrete Fourier transform (DFT) matrix 
and its conjugate transpose. Both CQVQ~^C and VQ~^CQV will be Hermitian and 
typically sparse, hence a number of dedicated eigenvalues solvers can be used to find 
the Slepian functions. In our own Matlab implementation we furthermore take ad- 
vantage of defining the operators as "anonymous" functions. This greatly improves 
the performance of the resulting algorithm, which we tested by comparison with the 
methods available for the special cases discussed in Sections |2] and [5] With this we 
have complete liberty to design Slepian functions on geographical domains, including 
the ability to construct "steerable", anisotropic windows for texture- sensitive analy- 



sis that is of interest in theory ( Olhede and Metikas^2QW\ \van de Ville and Unser 



2QQ8| ) and practice ( Audet and Mare schal 2007 ; Kirby and S wain 2006|). Fig. |7] shows 



two numerical examples of Slepian functions h spectrally sensitive in wedge-shaped 
oriented domains. The top row shows hi ^ h^ and their eigenvalues, Ai -^ A4, 
followed by eigenvalue- weighted sums of their periodograms, Xlc^^i Xa\Ha{'k)\'^, 
on a decibel scale. These are orientated at 7r/6 from the horizontal. The bottom 
row shows the equivalent results oriented at — 7r/6 from the horizontal. Both cases 
were computed by diagonalization (eigs) of the operator VQ~^CQV whereby Q is 
the anonymous function call to the two-dimensional Fast Fourier Transform (f f t2) 
and Q~^ its inverse (if ft 2). This immediately returns the /i(x) on a discretized 
calculation domain in which, once again, the region of the Colorado Plateaus was 
embedded. The functions hi ^ h^ shown are the real parts of the eigensolutions 
associated with the four eigenvalues Ai ^ A4 that have the largest real part. The 
eigenvalues of these "generalized prolate spheroidal data sequences" (gpss), which 
need not in principle be defined on regularly sampled grids (|^rwT£Z|1988j), are sen- 
sitive to the discretization and the size of the computational domains to which the 



inner products in eq. (91 ) refer. The slightly different manner in which these sam- 
ple the parametrically defined rotated spectral concentration domains explains the 
discrepancy of the eigenvalues between the two cases that are shown in the top and 
bottom rows of Fig. [7] 
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Fig. 7 Two sets of eigenfunctions /ii, /i2, . . . , /i4 that are strictly spacelimited within the Colorado 
Plateaus and optimally concentrated within a wedge-shaped spectral domain, with their concentration 
factors Ai, A2, . . . , A4. The last panel on each row shows the eigenvalue-weighted periodograms of the 
first q; = 1 — )► 20 eigenfunctions. The functions /i(x) = Vf{'x) were calculated by diagonalization 
of the spatial-domain operator equation {92) on a domain that we took to be three times the size of the 
box inscribing the concentration domain outlined in black, discretized on a 5 x 5 km^ grid. The spectral 
concentration domains are triangles with Hermitian symmetry oriented at ±7r/6 from the horizontal in the 
top and bottom rows, respectively. 



7 Conclusion 



Slepian functions are orthogonal families of functions that are all defined on a com- 
mon spatial domain, where they are either optimally concentrated or within which 
they are exactly limited. At the same time they are exactly confined within a cer- 
tain spectral interval, or maximally concentrated therein (Simons"20\Q). The joint 
optimization of quadratic spatio- spectral concentration is generally referred to as 
Slepian's problem, which we encountered as eq. ^ in one dimension, as eq. ( p3] ) 
on the surface of a sphere, and as eq. ([34]) in two Cartesian dimensions. Without 
qualification in one, and under special symmetry considerations in multiple dimen- 
sions, Slepian functions solve Sturm-Liouville equations, notably eqs. ( pTj ), ([24]), 
([26]) and ( [74] ). In those cases (though on the sphere only asymptotically) the solutions 
depend on the spatial and spectral areas of concentration only through the Shannon 
number, which depends on their product, as defined by eqs. ([9]), ( [2T] ) and ( [56] ). More 
generally, they solve Fredholm integral equations in their respective dimensions, as 
exemplified by eqs. ([5]), ( [TT] ), and ( [38] ). 

What unites the Slepian functions is that the kernels of the integral equations they 



solve are best thought of in the context of reproducing-kernel Hilbert spaces ( |Am/r- 
bekyan et al\2m^\N ashed and Walter\l99l\W^o\\96 1). Each of D(t, t') in eq. ([5bJ) 
D(r, y') in eq. \\lh) and I^(x, x^) in eq. ( [38b ), is a restriction of the traditional Dirac 
delta functions, ^(t, t'), ^(f , f '), and ^(x, x'), to a region of finite bandwidth {Simons 
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201Q| ). While the latter satisfy, for any square-integrable function /, the relations 

/oo 
f{t')S{t,t')dt', (93a) 

-CX) 

/(x)= / /(xO^(x,x0^x^ (93b) 

/(f) = / f{r^)5{r,r^)df2\ (93c) 

the former act on functions g that belong to appropriately bandlimited subspaces, as 
we have defined them in eqs. ([2]), ([14]) and ( [33] ), in much the same way: 



/CO 
g{t')D{t,t')dt', 
-CX) 

^(x)= / ^(xOI^(x,xo^x^ 

^(f)= f g{v')D{v,v')dQ. 



(94a) 



(94b) 



(94c) 



Note that these are not equal to eqs. ([5]), ( [TT] ), or ( [38] ). Hence the ^ are a basis 
for bandlimited processes anywhere on the applicable domain (R, R^ or the entire 
spherical surface i?) ( Dai/Z^^c/zI^ 1992[ |FtoJn>?||1998 Freeden et a/.||1998 Lan 
dau and Pollak\\96\\ Slepian and Pollak\1961 ). Therein, the Shannon-number best 



time- or space-concentrated members allow for sparse, approximate expansions of 
signals that are concentrated to the spatial region of choice. Similarly, the infinite 
sets of exactly time- or spacelimited (and thus band-concentrated) functions h which 
are the eigenf unctions of eqs. ([5]), ( [TT] ) and ( [38] ) with the domains appropriately re- 
stricted, see eq. ([45]) for the two-dimensional case, are complete bases for square- 



integrable scalar functions on the intervals to which they are confined { \Landau and 



Pollak\l96T\\Simons et al.\260S\\Slepian and Pollak\l96l} . Expansions of wideband 
but spectrally concentrated signals in the small subset of their Shannon-number most 
band-concentrated members provide approximations which are spectrally faithful and 
constructive as their numbers grow { Simons\20l0\ Simons et a/.|200 9). 

Slepian functions defined on two-dimensional Cartesian domains are of great util- 
ity in the natural sciences. Despite this, no comprehensive prior treatment was avail- 
able, beyond that of cases for regions with advanced symmetry. The present paper has 
attempted to do right by, especially, the geosciences, where Fourier-based methods 
and analyses are often required but the domains of interest are rarely circles or rectan- 
gles. Much can be learned from those special cases, however, including from the one- 
dimensional and spherical situations that we have reviewed also. Most particularly, 
they have helped us test and benchmark the algorithms that we have discussed and 
are making available on the Web as part of this publication. In having ended Sect. [6] 
with a computational procedure for the design of planar Cartesian Slepian functions 
on arbitrarily irregular spatial and spectral domains, which could be captured in just 



a few lines, see eqs. ([89])-([92]), and carried out in a handful of lines of Matlab code, 
we believe to have returned to a level of practicality that should appeal to researchers 
across a wide spectrum. 
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